The dynamics of large-scale structure formation in the 
universe by gravitational instability still presents many 
open issues and is mostly studied through numerical sim- 
ulations. This motivates the study of simpler models which 
can be investigated by analytical means in order to under- 
stand the main processes at work. Thus, we describe here 
a one-dimensional self-gravitating system derived from the 
cosmological context, which leads to an effective external 
potential in addition to the standard gravitational self- 
interaction. As a first step we consider small times so 
that the expansion can be neglected. Then we present a 
thermodynamical analysis of this system as well as the 
stability properties of the associated hydrodynamical and 

■ collisionless systems. We consider the mean field limit (i.e. 
continuum limit) to perform an analytical study. We find a 

■ second-order phase transition at T cl from an homogeneous 
| equilibrium at high temperature to a clustered phase (with 

Q-i a density peak at one of the boundaries of the system) at 
' low temperature. There also exists an infinite series of un- 

■ stable equilibria which appear at lower temperatures T cn , 
reflecting the scale-free nature of the gravitational inter- 
action and the usual Jeans instability. We find that, as 

\ for the similar HMF (Hamiltonian mean field) model, all 

■ three micro-canonical, canonical and grand-canonical en- 
' sembles agree with each other, as well as with the stabil- 

OS ity properties associated with a hydrodynamical approach. 
\ On the other hand, the collisionless dynamics governed by 
• the Vlasov equation yields the same results except that at 

\Q \ low T the equilibrium associated with two density peaks 

<^> ■ (one at each boundary) becomes stable. 
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Abstract. 



1. Introduction 

In standard cosmological scenarios the large scale struc- 
tures we observe in the present universe (such as clusters, 
filaments or voids) have formed through the amplifica- 
tion by gravitational instability of small primordial per- 
turbations, see Peebles (1980). Moreover, the amplitude 
of these density fluctuations increases at small scales as 
in the CDM model (Peebles 1982) which leads to a hi- 
erarchical process where smaller scales become non-linear 
first. Then, at large scales or at early times one can use 
a perturbative approach to study the evolution of ini- 
tial fluctuations (Fry 1984; Goroff et al. 1986; Bernardeau 
1992; Valageas 2001). Next, the weakly non-linear regime 
may be investigated through the Zeldovich approxima- 
tion (Shandarin & Zeldovich 1989) or the adhesion model 
(Gurbatov et al. 1989). However, the highly non-linear 
regime (which corresponds to collapsed structures such 
as clusters of galaxies) has mostly remained out of reach 
of systematic approaches. Thus, one may use the Press- 
Schechter approximation (Press & Schechter 1974) or the 
excursion-set formalism of Bond et al. (1991) to obtain 
the statistics of just-virialized halos or the saddle-point 
approach of Valageas (2002) for rare voids. These approx- 
imations focus on specific objects which may be "recog- 
nized" from individual features in the linear density field 
itself and do not follow the system as a whole through its 
non-linear evolution. A systematic method to do so was 
recently developed in Valageas (2004) but its application 
to the collisionlcss dynamics has not been performed yet 
for the highly non-linear regime. Therefore, the non-linear 
regime of cosmological structure formation is mostly stud- 
ied through numerical simulations. 

In order to simplify this difficult problem one can inves- 
tigate one-dimensional (1-D) systems which are easier to 
study both through numerical and analytical approaches. 
In fact, 1-D systems such as the system of parallel mass 
sheets (Camm 1950) have been studied for a long time to 
explore the evolution of isolated N— body systems which 
only interact through classical gravity. In particular, they 
have been used to investigate relaxation processes and to 



test Lynden-Bell's prediction for the final state of violent 
relaxation (Lynden-Bell 1967). Such 1-D systems can al- 
ready exhibit complex behaviours. For instance, Luwel 
& Severne (1985) found that collisional effects are not 
sufficient to relax the stationary waterbag configuration 
towards thcrmodynamical equilibrium whereas Rouet & 
Feix (1999) showed that holes in phase-space in the ini- 
tial distribution function can persist over long times and 
prevent efficient relaxation. On the other hand, the system 
can break up into smaller clusters as in Hohl & Feix (1967) 
or develop a fractal structure as in Koyama & Konishi 

(2001) . 

In this paper, we consider the 1-D gravitational sys- 
tem obtained by studying 1-D density fluctuations in a 3- 
D cosmological background. This model has already been 
studied mostly through numerical simulation in Aurell & 
Fanelli (2001), Aurell et al. (2001) and Fanelli & Aurell 

(2002) , both with and without cosmological expansion. 
Here we restrict ourselves to time-scales much smaller 
than the Hubble time so that the expansion of the uni- 
verse can be neglected and we focus on a mean field anal- 
ysis (i.e. a continuum limit) which is relevant in the cos- 
mological context. We study the thermodynamics and sta- 
bility of this system, similarly to the theoretical analysis 
of Chavanis et al. (2005) performed for the HMF model 
(defined by a cosine interaction) which shows a similar 
behaviour. We introduce this model in sect.El where the 
effect of the cosmological background is seen to reduce to 
an effective external potential V which balances the grav- 
itational self-interaction $ so that the homogeneous state 
is an equilibrium solution. Next, we perform a thermody- 
namical analysis of this system in sect. for the micro- 
canonical, canonical and grand-canonical ensembles. We 
study both the possible equilibrium states and their stabil- 
ity properties. Then, in sect. 0] we investigate the stability 
of such equilibria within an hydrodynamical framework 
whereas we consider the collisionless case as described by 
the Vlasov dynamics in sect. Finally, we present our 
conclusions in sect. El 
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2. A 1-D cosmological gravitational system 

The formation of large-scale structures in the universe 
through gravitational instability is described by the 
Vlasov-Poisson system which applies to a collisionless fluid 
in the continuum limit where the mass m of the particles 
goes to zero. Thus, the distribution function (phase-space 
density) /(x, v, t), where we note x the comoving coordi- 
nate of a particle and v its peculiar velocity defined by 
v = a 2 x (and a(t) is the cosmological scale-factor), obeys 
the Vlasov equation (see Peebles 1980): 



dt a 2 <9x <9x <9v 



(1) 



while the gravitational potential (j> satisfies the Poisson 
equation: 



—(p-p) with p(x,t = / /(x,v,t) dv. (2) 
a J 



Here we introduced the comoving density p(x, t) and its 
average p (which is independent of time) , and we note Q 
Newton's constant. The fact that the Poisson equation in- 
volves the density perturbation (p — p) rather than the 
total density p comes from the change from physical co- 
ordinates to comoving coordinates. In particular, for an 
homogeneous universe which merely follows the Hubble 
flow we have = and /(x, v,t) = p<5o(v) where Sd is 
Dirac's distribution. If we consider planar perturbations 
along the axis x: /(x, v, t) — f(x,v,t)5D(v±), we obtain 
a one-dimensional (1-D) problem: 



dt a 2 dx dx dv ' dx 2 



(p-p)- 



(3) 



Let us moreover restrict ourselves to periodic systems of 
period 2L which are symmetric with respect to the origin 
x = 0. Then, by symmetry the force F = —dfjj/dx van- 
ishes at the origin, F(x = 0) = 0, and the gravitational 
potential can be integrated as: 



2tt£ 



dx' (p — p){x) sign(x' — x), (4) 



and: 



dx'(p — p)(x ) \x — x'\, 



(5) 



where we restrict ourselves to < x < L. If we consider 
time-scales which are much smaller than the Hubble time 
H^ 1 = a/a we can neglect the time dependence of the 
scale- factor a(t) and taking a — - 1 and g = 2ttQ we can 
write the full potential 4> as: 



$ + V with 



$ = g J Q dx'p(x') \x - x'\ 



V(x) 



-gp 



(6) 



where we split <fi into a standard gravitational "interac- 
tion" term $ and an "external" potential V . Thus, the 



associated discrete system of N particles of finite mass m 
is described by the Hamiltonian Hjy: 



\ mv l + gm 2 \ x i - x . 

i=l 

N 

ftp/71 Y 



1>J 



I? 

4 



(7) 



Alternatively, the mean field approximation to the 
N— body system Q is given by the potential © and the 
Vlasov-Poisson equations. The important feature of this 
one-dimensional static cosmology (OSC) model which is 
due to its cosmological setting is the appearance of the 
"background" potential V(x) which ensures that the ho- 
mogeneous configuration p = ~p is a solution of the equa- 
tion of motion. This is the counterpart of the Hubble 
flow. The OSC model defined by JJJ or JJjJ was studied 
numerically in Aurell et al. (2001) and Aurell & Fanelli 
(2001) (where it is called a "frictionless model") for ran- 
dom (Brownian motion) and non-random (sinusoidal) ini- 
tial conditions. In this article, we shall consider the ther- 
modynamics of the continuous system 10, which can also 
be seen as the mean field approximation to (JJJ. This is a 
system of finite size L with < x < L. However, it can 
be useful to remember that it can also be understood as a 
periodic system of size 2L which is symmetric with respect 
to the origin x = 0. 

Here, we should add a few words about the limitations 
of this OSC model with regard to cosmology. The real 
universe is obviously three-dimensional and some impor- 
tant features such as the support of self-gravitating ha- 
los by rotation (as for spiral galaxies) cannot be studied 
through 1-D models. Thus, although the early stages of 
non-linear evolution are governed by 1-D collapse (forma- 
tion of pancakes, Shandarin & Zeldovich 1989) the highly 
non-linear regime which is of interest to us is truly three- 
dimensional. However, we can expect that key processes 
associated with the scale-free nature of the gravitational 
interaction are common to 1-D, 2-D and 3-D systems. For 
instance Munshi et al. (1999) find that 2-D and 3-D cosmo- 
logical simulations exhibit similar behaviours in the highly 
non-linear regime. In particular, 1-D systems should help 
to check whether the density profiles of virialized halos 
show a universal shape, for a wide class of initial condi- 
tions, and whether the power-spectrum of density fluc- 
tuations saturates into the non-linear regime or exhibits 
a power-law growth which could depend on initial con- 
ditions. A study of 1-D systems such as the OSC model 
is particularly timely as modern supercomputers are be- 
ginning to allow 1-D cosmological simulations in phase- 
space instead of the usual N— body techniques (Alard & 
Colombi 2005). 

On the other hand, the use of a finite size domain 
< x < L (which is equivalent to periodic bound- 
ary conditions of period 2L with symmetry with respect 
to x = 0) is customary to studies of gravitational sys- 
tems (e.g. Padmanabhan 1990). Within the cosmological 
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context, characterized by a hierarchical scenario where 
smaller scales turn non-linear first, the size L has actu- 
ally a physical meaning. It should be interpreted as the 
wavelength which turns non-linear at the time of inter- 
est. Indeed, larger scales which are still linear follow the 
Hubble flow (up to small perturbations) while smaller 
scales have relaxed within virialized halos. Thus, to this 
scale is also associated a characteristic time tn given 
by the age of the universe when L enters the non-linear 
regime. Today this yields tn ~ 10 10 yrs and L ~ 5 Mpc. 
Then, the expansion of the universe can be neglected if we 
restrict ourselves to length and time scales much smaller 
than L and tn- Of course, this may also prevent the sys- 
tem to reach thermodynamical equilibrium which can re- 
quire very long timescales (increasing algebraically with 
the number N of particles, Chavanis 2006). 

However, this article should only be seen as a first step 
in the study of the OSC system. We think that this model 
is also interesting in its own right as it should exhibit the 
features associated with a long-range scale-free interac- 
tion, which is an active domain of research. Besides, the 
background — p in Poisson eq.10 which arises from the 
cosmological context also leads to interesting transition 
phenomena as described in sect. below. 

3. Thermodynamical analysis 

3.1. Statistical equilibrium 
3.1.1. Micro-canonical ensemble 

Equilibrium profiles: 

In the micro-canonical ensemble the statistical equi- 
librium is obtained by maximizing the Boltzmann en- 
tropy S[f] at fixed energy E and mass M, given by 
(Padmanabhan 1990, Chavanis et al. 2005): 

S = -Jdxdvf\nf, E = Jdxdvf(j + ~+v) (8) 



and: 

M = pL= I dxdvf. 



(9) 



Using Lagrange multipliers this yields the usual Maxwell- 
Boltzmann distribution: 

/(i,!))oce-* 2 / 2+ ^ and p(x) oc e ~W*\ (10) 

where we introduced the inverse temperature = 1/T. 
Let us define the offset Tp such that: 



ip + <p with p — pe P v and — j — 2g(p — p). (11) 



dV 

dx 2 



Therefore, the thermodynamical equilibrium is set by the 
equations: 



ay 

dx 2 



where the boundary conditions are obtained from eq.Q. 
Without the constant term —1 in its r.h.s. ca. (|12ll re- 
duces to the usual Lane-Emden equation which arises for 
isolated self-gravitating isothermal gaseous spheres (e.g., 
Chandrasekhar 1942; Binney & Tremaine 1987). It is 
convenient to introduce the Jeans length Lj (Binney & 
Tremaine 1987) and the ratio (l- 

Lj = (2gppy 1/2 , & = ^ = y/^L, (13) 



in order to define the dimensionless quantities: 



C = yd,, 4> = P<P and rj = = , 
T, p 



— zi 7 whence r\ 
P 

Then, the modified Lane-Emden ea. I12|) reads 



d 2 ip 
dC 2 



1 and V'(0) = iI>'(Cl) = 0. 



= e - * (14) 



(15) 



This is the equation of motion of a particle of mass m = 1 
and coordinate ip as a function of time £ within a potential 
V(ip) — + ip — 1, with zero velocity at initial (£ = 0) 
and final (£l) times. The case where the particle is at rest 
at the bottom ip = of the potential V(ip) is always a so- 
lution of eq. (|15|l . This corresponds to the uniform density 
p = p. As the temperature T decreases, while (3 and (l 
grow, new solutions appear as the particle now has time 
to perform one or several oscillations in the potential well 
V(ip). This yields a series of critical temperatures T cn , or 
(3 cn , associated with the appearance of a solution with n 
half-oscillations. They can be obtained by expanding the 
potential V(ip) around ip = up to quadratic order. This 
yields: 

ib 2 

V(V0 — — j whence ip" + tp = and ip — ip cos(, (16) 

where we used ip'{Q) = 0. Next, the second constraint 
V>'(Cl) = implies: 



Cl = nir whence j3 cn = 



2gpL 2 



T — 



2gpL 2 



(17) 



2gp (e~ 0v ~ 1) and <p'(0) = tp'(L) = 0, (12) 



If n is odd this corresponds to n half-oscillations starting 
at tpo and ending at —ipo whereas if n is even it corre- 
sponds to n/2 full oscillations starting and ending at ipo- 
As the temperature decreases below T cn the oscillations 
get wider and take "more time" . In the following we shall 
label by "n" the solution with n half-oscillations which 
starts at ^(0) = ip- < and by — n the solution with 
n half-oscillations which starts at ^>(0) = ip+ > 0. From 
the last equation in (|10fl or (|14l) we see for instance that 
the mode n — 2 shows two density peaks at x — and 
x — L while the mode n — —2 shows only one density 
peak at x — L/2 (but two underdense minima at x = 
and x = L). We shall call the homogeneous equilibrium so- 
lution p = ~p as the mode n = 0. On the other hand, since 
the differential equation in l|15|) does not depend on fj or n 
(which only appear in the boundary conditions) all modes 
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|8 = 4j8 al 
/8 = 2/Soi 



Fig. 1. The equilibrium profile n = 1 obtained at /? = 2/3 c i 
(dotted lines) and /3 = 4/3 c i (dashed lines). The solid lines 
correspond to the uniform solution (n — 0) with rj = 1 and 
•0 = 0. The normalized density profiles r/ are the curves 
with a maximum at £ = while the normalized potentials 
V> have a maximum at £ = (i. 

±n can be written in terms of the mode n = 1 at a dif- 
ferent temperature Ti given by £li = Ct/ n or A = /3/n 2 . 
Thus, we have: 



MC;P) = ix C 



V-„(C;/?) = (C-— ;4 



(18) 
(19) 



and similar relations for the overdensity rj. Here we con- 
sidered that all solutions are extended over the full real 
axis to yield solutions which are periodic of period 2L and 
are symmetric with respect to x — 0, in agreement with 
sect. |2 This merely amounts to follow the motion of the 
particle within the potential well V(ip) over all times. The 
periodicity and symmetry of "0(C) implies that we can ex- 
pand the solution as: 



i(0 



2 O-k cos I KIT— 

fe=0 ^ ^ L 



(20) 



This automatically satisfies the boundary conditions l|15|) ■ 
As seen in eas. (ll8[l - H19|) it is sufficient to consider the case 
n = 1. Then, in order to study the behaviour close to the 
critical temperature we can substitute the expansion (1201) 
into the differential ea. 1)15(1 up to the required order over 
{(3 — (3 c i). We obtain up to order (/3 — /3 c i): 



a-2 



12' 



a i 



'12 



- Pc 



hi 



(21) 



for k > 1. 



while higher-order terms are of order ' 
Thus we obtain a square-root singularity for the amplitude 
of the potential ip and for the density profile. 



On the other hand, we can also obtain asymptotic ex- 
pressions for the equilibrium density profile in the limit 
of low temperatures T — ► 0. Thus, close to ( = where 
ip — > — oo we can approximate V(ip) — which yields 
from ea. ((15(1 



InCr 



c 



which leads to 

In a 



dip' 



Ci 



v(0 



Cl 



2cosh' ! (CLC/2) 



(22) 



(23) 



As expected, we recover the equilibrium profile obtained 
by Camm (1950) for an isolated 1-D gravitational sys- 
tem. Indeed, close to the density peak we can neglect p in 
Poisson equation ((11(1 in the low temperature limit (since 
p(0) r~j 1/T from ea. ((23(0 . Thus for arbitrary n we obtain 
a series of density peaks which follow the equilibrium pro- 
file of isolated 1-D gravitational systems separated by low 
density regions since in the case of the mode n = 1 we 
obtain: 



3 -tt 2 T c i/2T 



(24) 



where we used the approximation V(ip) — ip in the un- 
derdense region. We can note that the density in these 
underdense regions decreases very fast at low tempera- 
tures. 

We display in fig. ^ the equilibrium profiles n = 
(homogeneous solution) and n = 1 (one density peak at 
C = 0) at temperatures T c i/2 and T c \/&. We can already 
check at these temperatures that the density peak grows 
as 1/T while its width decreases as T and that the un- 
derdense region is almost void. Of course, the potential 
ip = — Inr) shows a much shallower shape. 

Thermodynamical quantities: 

From ea. ((10(l the distribution function at equilibrium 

is: 



f(x,v)=p(x)J^ e~^l\ 



(25) 



whence the entropy S defined in eq.JHJ is up to an addi- 
tional constant: 



M 



InT- 



da; pin ^ = — InT + p I dxr/ip, (26) 
p 2 



where we used eas. (11411 . Therefore, the specific entropy s 
normalized such that s(T c i) = is: 



S - S(T c i) 1 T 

s = = — In — 

M 2 T cl 



Cl 



rjip. 



(27) 



On the other hand, the energy E at thermodynamical 
equilibrium can be written from eq.JSJ and ea. ((10(l as: 



E = 



MT 



d*p— with p = p+- — . 



(28) 



6 



P. Valageas: Thermodynamics and dynamics of a 1-D gravitational system 





Fig. 2. The specific entropy s = (S - S{T cX ))/M as a 
function of the energy E for the modes n = 0, ±1 and ±2. 
The lower energy bound E ni i n (n) associated with mode 
±n, reached at T = 0, increases for higher n > 0. The ho- 
mogeneous solution n — has the highest bound En. All 
thermodynamical quantities are identical for both equilib- 
ria n and —n. 



Fig. 3. The specific heat C v — dE/dT as a function of 
temperature for the equilibria n — 0, ±1 and ±2. It is 
discontinuous at the transitions T cn and corresponds to a 
second-order phase transition. 



Integrating over the term pV, using the second eci. (|28J) 
and integration by parts, yields: 

MT gM 2 L M, , N , „ /" , , , 

-^^ + T [0(O) + 0(i)] + y dx (p-p)|.(29) 

From eq.JSJ) we obtain </)(0)+</>(i) = therefore the energy 
of modes and ±n reads: 



To 



MT gM 2 L 
~2 ~ 



(30) 



MT gM 2 L MT [ a dC , 



6 



Fig. 4. The calorific curve /3(E) for the modes n = 0, ±1 
and ±2. It shows a second-order phase transition at T c i. 



In particular, the energy E cn associated with the critical 
temperature T cn of ea. 1)17(1 is: 



jr. — 



n it — 6 
6n 2 7r 2 



gM 2 L at T c; 



2gML 



(32) 



Here we can note that the energy E of the OSC system 
is actually bounded from below. Indeed, as can be clearly 
seen from eq.(|7|) the minimum value E m - ln for the energy 
is reached when all the mass is at rest at x = or x = L. 
On the other hand the energy En associated with the 
homogeneous static solution p = ~p and v = 0, which would 
be the counterpart of the original Hubble flow, can be 
obtained from ea. (|3(J|l with T = 0. This yields: 



gM 2 L 



and Th = — 



gM 2 L 

6 



(33) 



We can note that since Th = Eq(T = 0) all equilib- 
rium solutions ±n can be reached with this energy, as 
may also be seen from fig. [21 where we display the specific 
entropy s as a function of the energy E. As we can check in 
fig. [21 the mode n = ±1 reaches the energy bound E min at 
T = while other equilibrium states ±n are increasingly 
close for larger n to the homogeneous solution, with re- 
spect to global quantities such as energy or entropy, since 
they are made of n/2 density peaks separated by a length 
2£i/n, and their minimal energy -E m i n (n) gets closer to 
En. Although fig. [21 only shows the modes 0,±1 and ±2 
we can see that at a fixed energy E the homogeneous 
state has the lowest entropy while modes ±n with n > 1 
have a higher entropy as n decreases. In particular, below 
T cl the mode ±1 with only one density peak at x = 
or x = L is the highest entropy equilibrium. Therefore, 
it is a stable equilibrium solution of the system, as we 
shall also check in sect. 13.2.11 whereas the homogeneous 
solution is unstable. Thus, as in the case of isolated 1-D 
systems (where only equilibrium ±1 exists and extends to 
all temperatures) and contrary to 3-D systems we have 
no gravothermal catastrophe (associated with the disap- 
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pearance of any local entropy maximum, e.g. Binney & 
Tremaine 1987; Padmanabhan 1990). 

Of course, the energy E is continuous at the transitions 
T cn . However, its derivative C v = AE/AT is discontinuous 
as we jump onto the branch ±n. Indeed, from ea. H30|l the 
specific heat C„o associated with the uniform solution is: 



M 



(34) 



while near the transition T cl , using the expansion Q2U[1- 
J2U, we obtain: 

7M 7M 
E x ~E cl + —(T-T cl ) hence <7 i(I£) = — ■ (35) 

We show in fig. the specific heat for n = 0, ±1 and ±2. 
We can check the discontinuity at T cn which corresponds 
to a second-order phase transitions. Next, we display in 
fig. 0] the calorific curve /3(E). It also clearly exhibits this 
second-order phase transition and we can see on the fig- 
ure the various low energy bounds E m i n (n). Besides, we 
can see in fig. 0] that the specific heat is always positive 
and no gravothermal catastrophe can develop. Overall, 
the behaviour of the OSC model is quite similar to the 
HMF model studied in details in Chavanis et al. (2005). 
The latter model is defined by a cosine self-interaction 
cos[27r(:Ei — Xj)/L] instead of \xi — Xj\ without the "ex- 
ternal" potential V. It can also be seen as a truncation 
of the 1-D gravitational potential $ to its first Fourier 
component. As described in Chavanis et al. (2005) the 
HMF model also exhibits a second-order phase transition 
from an homogeneous equilibrium above T cl to a clustered 
phase with only one density peak below T c \. However, 
there are no other (unstable) equilibria ±n appearing at 
T cn . Of course, the infinite series of equilibria ±71 which 
appear in the gravitational case is due to the scale-free na- 
ture of the gravitational potential cx \xi— Xj\ (which is only 
broken by the finite length L of the system which leads to 
a discrete spectrum instead of a continuous one). We shall 
check in the following that apart for the additional equi- 
libria ±n the properties of the HMF system and of the 
OSC model are similar both from thermodynamical 
and dynamical points of view. 

3.1.2. Canonical ensemble 

In the canonical ensemble the thermodynamical equilib- 
rium is obtained by minimizing the free energy F = 
E — TS at fixed temperature T and mass M. This yields 
again the Maxwell-Boltzmann distribution (|10|) hence we 
recover the same critical temperatures and solutions (|17|) - 
The free energy F(T) reads: 



F + TS(T cl 
MT 



In I — 

2 \T cl 



E-MTs 
MT 



MT gM 2 L 



2 

Cl AC 



6 

(f7+l)V. 



(36) 




Fig. 5. The free energy f(T) = F(T) + TS(T cl ) as a func- 
tion of temperature for the states n = 0, ±1 and ±2. 

We can see that at fixed T the homogeneous state (n = 0) 
has the highest free energy while F decreases for lower n 
with n > 1. Therefore, as we shall check in sect. 13 .2. "21 the 
uniform equilibrium is again unstable below T c % while the 
equilibrium n = ±1 is stable. 

3.1.3. Grand-canonical ensemble 

In the grand-canonical ensemble the thermodynamical 
equilibrium is obtained by minimizing the grand potential 
fi = E — TS — fiM at fixed temperature T and chemical 
potential p. Since the mass M now fluctuates the mean 
density (p) = M/L within [0,L] is no longer set equal to 
p. In this case, we cannot any longer extend the system to 
the whole real line using a periodicity of 2L and symme- 
try with respect to x = 0. It can only be interpreted as a 
small system of actual size L within a larger cosmological 
background of mean density p, analyzed over time-scales 
much smaller than the Hubble time. Thus, if (p) > ~p it 
corresponds to a local overdensity while if (p) < ~p it cor- 
responds to a local underdensity. Then, minimizing the 
grand potential fl yields the same results (|17J) - H21J) as for 
the micro-canonical and canonical ensembles. Note that 
the agreement between all three thermodynamical ensem- 
bles does not hold for all 1-D systems. For instance a sys- 
tem of concentric spherical cells which interact by gravity 
can show metastable states, negative specific heat and dis- 
crepancies between different thermodynamical ensembles 
as seen in Youngkins & Miller (2000). Thus the 1-D sys- 
tem studied here from planar perturbations shows a rather 
simple behaviour while presenting a phase-transition and 
a scale-free origin. 

3.2. Thermodynamical stability 
3.2.1. Micro-canonical ensemble 



We display in fig.0the free energy F as F(T)+TS(T cl ) as 
a function of temperature for equilibria n = 0, ±1 and ±2. 



As recalled in sect. !3.1d1 the thermodynamical equilibrium 
states within the micro-canonical ensemble are obtained 
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by maximizing the entropy S. In sect. 13.1.11 we investi- 
gated the first variations of the entropy S, which yields 
maxima, minima and saddle-points. Next, the thermody- 
namical stability of these states is obtained by studying 
the second variations of the entropy. A stable equilib- 
rium state corresponds to negative definite second vari- 
ations. Following Padmanabhan (1990), maximizing the 
entropy S[f] at fixed energy E and density p(x) yields the 
Maxwellian distribution Q25JI. Then, we can write the en- 
tropy S and the energy E in terms of the density p(x) as 
in eas. (|26|) . (|28[) . Next, we can obtain from these expres- 
sions the second variation S 2 S with respect to a density 
fluctuation 5p{x) at fixed energy (the constraint SE = 
yields 8T in terms of 5p, which is next substituted into 
SS). This yields for the second variation around an equi- 
librium solution: 



a2 c 9_ 

6 b ~ 2T 



MT 2 



dxdx' 5p{x)8p(x')\x — x'\ 
dx <f) 8p 



dx 



6 Z 

2p 



(37) 



As in Padmanabhan (1990) let us define the mass fluctu- 
ation q(x) below x as: 



OP = 

dx 



q(x) 



where we used the conservation of mass M. Then, after 
integrations by parts the second variation S 2 S can be writ- 
ten as: 



S 2 S = / dxdx' q(x) Kmc {x,x')q{x') 
Jo 



(39) 



with: 

K M c{x,x') = —5 D {x-x')- 



T 



MT 2 dx 1 'dx { ' 



ls D (x-x') — 



1 



dx' \p(x') dx' 



(40) 




Fig. 6. The largest two stability eigenvalues Xk — Xk/g/3 
(i.e. = 1,2) as a function of inverse temperature (3 
for the equilibrium states (solid lines), ±1 (dashed 
lines), 2 (dotted lines) and —2 (dot-dashed lines) for the 
micro-canonical ensemble. Stability corresponds to nega- 
tive eigenvalues Xk < 0. The states ±1 are actually iden- 
tical (up to a symmetry with respect to x = L/2). The 
equilibria ±n show n — 1 unstable eigenvalues while the 
homogeneous solution n — shows n unstable eigenvalues 



dx' Sp{x'), g(0) = q(L) = 0, (38) below T c , 



Thus, the equilibrium p(x) is stable if all eigenvalues A of 
the kernel Kmc are negative (so that S[p] is a maximum). 
Hence we are led to investigate the eigenvalue problem 
Kmc-1 = Ag which reads: 



l__d_ 

2dx 



1 dq 



p(x) dx 



9 ( \ 



Xq(x) 



MT 2 



dx'(j>'q. (41) 



For the homogeneous equilibrium (n = 0) where ?y = 1 
and '0 = we obtain: 

q" + (1 - X)q = 0, g(0) = g(a) = 0, (44) 
whose solutions are: 

q k = A k sin (kn^j , X k = 1 - (^j , k = 1, 2, .. (45) 

Thus, as expected we find that the uniform equilibrium is 
stable at high temperature T > T c \ (i.e. £l < 7r). At each 
critical temperature T cn defined in eci. (|17f) a new mode 
of instability A„ appears while two new equilibrium states 
±n appear which exhibit n— 1 unstable modes. This agrees 
with fig.[21which displays the specific entropy as a function 
of energy. In particular, close to T~ x we can investigate the 
stability of equilibrium n = ±1 by using the expansion 
H20fl - H21fl and looking for a similar expansion for X\. Thus, 
we write: 



In terms of the dimensionless variables (|14fl we obtain the 
eigenvalue problem: 



a? 9i,2 



Ai = a 1 Ai i 



(46) 



_d_ 



1 dq 



q = Xq + 2i>' I ^'q, with A = A, (42) 

Cl g/3 



which yields (we can check that Ai.i = 0): 



supplemented by the boundary conditions q(0) = g(Ci) = 
0. We can note that the matrix Kmc defined by: 



K M c-q = - A 



1 dq 

vd( 



■><•' I ^'q 



(43) 



is positive definite, therefore A < 1. 



1 5 
9i.o = sinC, 91.1 = B x sinC -- sin2C, Ai, 2 = -j- (47) 

Therefore, since Ai < close to T~ x we see that equilibria 
±1 are stable close to T~ x while equilibrium becomes un- 
stable. From fig. |21 we see that the states ±1 are actually 
stable down to T = 0. We display in fig.[|)]the first two sta- 
bility eigenvalues Xk as a function of /3 for the equilibrium 
solutions n — 0, ±1 and ±2, obtained from a numerical 
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Fig. 7. The stability eigenvalues = \k/gfl as a function 
of inverse temperature (3 for the equilibria 0, ±1 and ±2 for 
the canonical ensemble. Stability corresponds to negative 
eigenvalues A& < 0. 



Fig. 8. The stability eigenvalues Afc = Xk/g/3 as a function 
of inverse temperature (3 for the equilibria 0, ±1 and ±2 
for the grand-canonical ensemble. Stability corresponds to 
negative eigenvalues A& < 0. 



study of eq.|@2J. We can check that the equilibria ±1 are 
stable above (3 C \ while the homogeneous equilibrium is un- 
stable. Higher states ±n exhibit n— 1 unstable eigenvalues 
therefore they are only saddle-points of the entropy S[f] 
and they are not thermodynamically stable. 

The series of equilibrium states ±n which appear as 
new unstable modes from the homogeneous solution at 
increasingly low T cn are somewhat analogous to the high 
order modes of instability which appear for a 3-D finite 
isothermal sphere obtained in Semelin et al. (1999) and 
Chavanis (2002). They can be related to the scale- free na- 
ture of gravitational interactions (although the homoge- 
neous background yields an arithmetic hierarchy of scales 
while the isothermal sphere leads to a geometric progres- 
sion). However, since the equilibria ±n with n > 2 are un- 
stable (except n = 2 at low T for a collisionless dynamics 
as seen in sect.EJbelow) a discrete N— body system which 
undergoes collisional relaxation is not expected to frag- 
ment into many clumps at low T but to form only one col- 
lapsed density peak close to one boundary (though several 
temporary clumps might appear in a transient regime). 

3.2.2. Canonical ensemble 

Within the canonical ensemble we need to minimize the 
free energy F. Hence stable equilibria are characterized 
by negative definite second variations of — F at fixed tem- 
perature and mass. Following the same procedure as in 
sect. I3.2~T1 we obtain a similar eigenvalue problem with a 
kernel Kc(x,x') which is equal to Kmc{x,x') written in 
ea. l|injl without the term 4>'(x)(/)'(x')/MT 2 . Therefore, the 
eigenvalue problem now reads: 



d 



1 dq 



+ q = \q, q(0) = q(( L ) = 0. 



(48) 



We again obtain as in eq. ||4H) that A < 1. Then, the stabil- 
ity properties of the uniform solution are exactly the same 



as for the micro-canonical ensemble (since '0 = 0). On the 
other hand, for equilibrium 1 we can again use the expan- 
sion close to T~ v This now yields Ai,2 = -1/4. 
Therefore, states ±1 are again stable close to T~ x . We dis- 
play in fig. [3 the stability eigenvalues obtained from a 
numerical study of eq.Q48)l. We obtain a behaviour which 
is similar to fig. and we recover the same thermody- 
namical stability properties. This agrees with fig. [5] for 
the free energy. Note that the equivalence between the 
micro-canonical and the canonical ensemble (with, how- 
ever, different eigenvalues) obtained in this case is similar 
to that found for the HMF model but differs from the case 
of 3-D self-gravitating systems where the limits of stabil- 
ity do not coincide in micro-canonical and canonical en- 
sembles (Padmanabhan 1990; Chavanis 2002). This is also 
related to the fact that the specific heat is always positive 
(see fig. ISJ) which is not the case for 3-D self-gravitating 
systems (in the micro-canonical ensemble). 

3.2.3. Grand-canonical ensemble 

Within the grand-canonical ensemble we now need to min- 
imize the grand potential fi. The second variation of fl 
reads: 

„2 



J dxdx' Sp(x)S P (x')\x - x'\ +T J dx (49) 

Since the total mass is now allowed to fluctuate we define 
the perturbation 8~p of the mean density M/L and q(x) 
by: 

Sp(x)=Sp + q'(x), q(0) = q(L) = 0. (50) 

Substituting into ea. (|49J) and integrating by parts yields: 

r T 1 r a / V \ £— 2 

s 2 n = -q I dx q 2 -- 



dxq dx-{^) + — — 



dx q 



-T^ + g(2x - L) 
P 



(51) 
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Eq. l|51|) shows that 5 2 tt has a finite minimum over Sp. 
Thus, we first minimize <5 2 over Sp and the second vari- 
ation of with respect to q(x) alone now reads: 



s 2 n 1 
-— (?) = 2 



9jt 
3T 



dx 
P 



dec q 



.(52) 



This can again be written in terms of a kernel Kqc( x : x ') 
which leads to the eigenvalue problem: 



However, this analysis does not yield the growth rates as- 
sociated with such instabilities nor the dynamical proper- 
ties of the system. Within a mean field or continuum ap- 
proach we shall study in sect. El the dynamics as given by 
the Vlasov-Poisson system. However, before we tackle the 
Vlasov equation in phase-space (x,v) we first investigate 
in this section the hydrodynamical model associated with 
the hamiltonian Q . Indeed, the fluid and collisionless sys- 
tems often share important properties while the hydrody- 
namical system is much easier to analyze (e.g. Binney & 
Tremaine 1987). 



1 dq 
77 dC 

dc U 



(C-|) 



(C 



dC 
v 



y)) 



(53) 



For the uniform solution this equation reads: 



q" + (1 - X)q 



2 



dC(C'-^)tf(.'). <--n 



For odd k we recover the modes <|45[) of the micro- 
canonical and canonical ensembles (since the r.h.s. in 
ea. <|54|) vanishes) whereas even modes are modified. In par- 
ticular, their eigenvalues are now given by: 



8 



20 

with A = I - ( — ) .(55) 



For large 9 this yields 62k — kir hence we recover the 
even modes \2k — 1 — (2fc7r/^) 2 of l|45|). One can see 
from ea. H55l) that A is increased at fixed C,l as compared 
with the micro-canonical and canonical result l|45[l . We 
can note that for Ql > ^60 + 12V30 the first even mode 
has A2 > 1, that is 62 is imaginary and can be written as 
2 = i§ 2 with <9 2 coth(9 2 = 1 + el/2, - <9|(16/Cl + 8 / 3 C£)- 



In any case, we find that as for the micro-canonical and 
canonical ensembles the homogeneous equilibrium is un- 
stable below T c i while other equilibrium solutions ±n are 
also unstable except for n = 1 , as can be seen from fig. [S] 
where we show the results of a numerical analysis of the 
eigenvalue problem (J53J. Again, it is possible to analyze 
the stability of state ±1 near the transition T c \ using the 
expansion pU ]) -l|2T )l . This yields Aj.,2 = -1/4 < which 
implies that equilibria ±1 are again stable close to T~ x . 
From fig.|Hlwe see they are actually stable down to T = 0. 
Thus, all three thermodynamical ensembles lead to the 
same stability properties. 

4. Hydrodynamical model 

The thermodynamical analysis of sect. |3]has allowed us 
to obtain the equilibrium states of the OSC system (0 
and to derive their thermodynamical stability properties. 



4.1. Hydrostatic equilibrium 



The equations of motion of the hydrodynamical system 
associated with the OSC model JjJ are simply: 

dp . d 

at 



+ —(pv) - 0, 



dx 



(56) 



dv 

at 



dv 
dx 



1 dP 

p dx 



dx ' 



(57) 



where the potential <j> is given by eq.JHJ and P{p) is the 
barotropic pressure of the gas. As usual, one can check 
that the energy functional W is conserved by the equations 
of motion (Binney & Tremaine 1987), with: 



W = dx p(x) 



P 



dx p- 



+ + V 



(58) 



Therefore, minima of W are stationary solutions of the 
Euler equations which are non-linearly stable. As in the 
statistical analysis of sect. 13.11 the first variation of W 
yields the properties of such stationary states (in addition 
to saddle-points of W). Then, the first variation of W with 
respect to v leads to v = 0. Substituting into the Euler 
ea. (|57|) yields the hydrostatic equilibrium condition: 



(59) 



dP _ dej) 
dx ^ dx 

In the following we shall restrict ourselves to the case 
of the isothermal gas defined by the equation of state 
P = pT. Then ea. (|5^|l reads dp/dx = —flpd4>/dx which 
yields again p oc e - ^ as in the thermodynamical approach 
(|10|l . Therefore, we recover the same stationary states as 
in sect. 13.11 

4.2. Non-linear stability 

The equilibrium states obtained in sect. 14.11 are non- 
linearly stable if they are minima of the energy functional 
W, that is if the second variations of W are definite posi- 
tive. The second variation with respect to v is clearly pos- 
itive from the expression (|58|l while 5 2 W/Sp5v = since 
v = for equilibrium states. Therefore, we only need to 
study the second variation with respect to the density. For 
the isothermal gas this yields: 



S 2 W = 



1 



dec 5p5$ 



dx 



dp 2 



(60) 
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Fig. 9. The stability eigenvalues Uf. = \\/(2g~p) as a func- 
tion of inverse temperature (3 for the equilibria 0,±1 and 
±2 for the hydrodynamical system. Stability corresponds 
to imaginary eigenvalues \\ < 0. The second eigenvalues 
A2 obtained for equilibria ±2 are equal. 



This is the same equation as the second variation of the 
free energy F obtained for the canonical ensemble in 
sect . 13. P| Hence the non-linear stability properties of var- 
ious equilibrium states are the same as those derived for 
the canonical ensemble. Therefore, the non-linear hydro- 
dynamical stability properties are identical to the thermo- 
dynamical ones. This implies in particular from ea . (|45|l 
that the homogeneous equilibrium is unstable below T c \ 
while equilibria ±n are unstable close to T cn for n > 2. 
From fig. they remain unstable down to T — ► while 
equilibrium ±1 is stable below T&. 

4.3. Dynamical linear stability 

The analysis of sect. 14.21 allowed us to discriminate be- 
tween stable and unstable states but it does not provide 
the growth rates of unstable perturbations. To this order 
we need to investigate the dynamical linear stability of 
such equilibrium states. Thus, we linearize the equations 
of motion l|5fi [l -l|57 [l and we study the normal modes of 



the form Sp < 



Introducing again the mass fluctuation 



q(x) as in ea. (|38|) and using the dimensionless variables 
(|14|) we obtain the eigenvalue problem: 



d_ 



1 dq 



v— with v 
V 



2gp 



and Sv 



-Al(61) 
P 



Thus, eigenmodes with v < (i.e. A is imaginary) are sta- 
ble while eigenmodes with v > (i.e. A is real) are unsta- 
ble. More precisely, for v > we have both a decaying and 
a growing mode (A = ±\Jv2g~p) distinguished by opposite 
velocity fields, as seen from the last equation in (|61fl . At 
the point of marginal stability v = we recover the eigen- 
value problem l|48|) of the canonical ensemble. Therefore, 
from sect. l4.2l we see that the conditions of non- linear and 
linear stability are identical and also coincide with the 



Fig. 10. The unstable eigenmode k = 1 with v\ ~ 1 for 
the equilibrium state n = — 2 at temperature /3 = 6/3 c i. 
We plot the equilibrium density profile ij = p/~p (dashed 
line), the unstable eigenmode Srj = Sp/p (solid line) and 
the velocity Sv oc — q/rj (dotted line). This eigenmode cor- 
responds to a mere translation of the density peak. The 
normalizations of Srj and Sv are arbitrary. 



canonical ensemble. Note that this equivalence also holds 
for the HMF model (Chavanis et al. 2005) as well as for 
some more general models (Chavanis 2006). 

For the uniform equilibrium (n = 0) we have 7/ = 1 so 
that we recover the eigenvalue problem l|44(l of both the 
micro-canonical and canonical ensembles with the same 
solutions (|45(l . Thus we find again that the homogeneous 
solution is unstable below T c \. For the equilibrium state 
±1 we can again use the expansion 

J20J-I211) close to T cl 
which yields v\ ~ — af/4 < 0. Therefore the equilibrium 
±1 is again stable close to T c y. From fig. El which shows 
the eigenvalues Vk as a function of inverse temperature (3 
we can check that the equilibrium ±1 is stable below T c \ 
down to T = whereas equilibria ±n with n > 2 are al- 
ways unstable. This is identical to the non-linear stability 
properties obtained in sect. 14.21 and to the thermodynami- 
cal results of sect. 13.21 We can note that the equilibria ±2 
have different first eigenvalues v\ > but identical second 
eigenvalues v-i < 0. Indeed, we can check numerically that 
this eigenmode vanishes both at the boundaries x = 0, L, 
and at the midpoint L/2 (q is antisymmetric with respect 
to L/2) so that both rj and Srj are identical for n = 2 and 
n = — 2 if we extend them to the whole real axis up to 
a mere translation of L/2. This is due to the symmetry 
properties of equilibrium profiles ±2. 

Next, we can note that the eigenvalue v\ obtained for 
the state n — — 2 quickly becomes very close to unity 
above T c i- In fact, we can check that q = r\ is a solu- 
tion of the differential eq. l|t)ll) with v — 1, using r\ = 
and ip" = r; — 1 which characterize the equilibrium pro- 
files, see eas. lfT4^l and (|TK)l . Because of the boundary con- 
ditions q(0) = q(Cl) = this is not an allowed solution 
to the eigenmode problem l|tjlfl . However, for the equilib- 
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Fig. 11. The unstable eigenmode k = 1 with Vy > for 
the equilibrium state n = 2 at temperature (3 = 6/3 c i . We 
plot the equilibrium profile r\ (dashed line), the unstable 
eigenmode 8r\ (solid line) and the velocity 8v (dotted line) . 



rium n = — 2 where the overdensity goes to zero at the 
boundaries as rj ~ e~ v Tc2 / 2T from eq.||22J the unstable 
eigenmode tends to q oc rj at low T. From the definition 
this implies for the density perturbation Srj oc r/(C) 
which corresponds to a mere translation of the equilibrium 
density profile. We display in fig. the first eigenmode 
v\ ~ 1 for the equilibrium state n = — 2 at temperature 
T C 2/6. We plot the equilibrium profile 77 (dashed line) and 
the growing mode Srj (solid line) with the associated ve- 
locity Sv (dotted line). The derivative rj' (£) normalized 
to the Srj at its peak cannot be distinguished from the 
perturbation Srj on the figure and we can see that the ve- 
locity is almost constant over the extent of the density 
peak. Thus the unstable mode obtained for the equilib- 
rium n = — 2 only corresponds to a translation of the 
density peak (which behaves as an isolated system as dis- 
cussed in sect. 13.1.10 while other eigenmodes k > 2 are 
decaying modes, as seen from fig. El The associated eigen- 
value v\ — 1 merely corresponds to a point-particle which 
rolls down the "external" potential V(x) of eq.©. Indeed, 

the solution of x = — dV/dx reads (x — L/2) oc e 



±j2gpt 



which describes the motion in the low-temperature limit 
T — > as the density peak (|23|) tends to a Dirac. It would 
be interesting to follow the dynamics further from x = L/2 
to see whether the reflection at the boundaries x = 0, L, 
quickly destroys the equilibrium profiles or whether some 
stationary cycle can be reached (if there were true solitonic 
solutions). We shall leave this analysis for future work as 
it goes beyond the scope of the present paper. 

Finally, we can see in fig. [5] that for the equilibrium 
n — 2 the eigenvalue vi quickly goes to zero. Of course, 
we cannot recover the solution q oc 77 associated with a 
mere translation since the state n = 2 is made of two den- 
sity peaks at x = 0,L, which grow at low T so that the 
boundary conditions q(0) = q(L) = forbid the solution 
q oc 77. We show in fig. ^2 the first eigenmode v\ > for 



the equilibrium state n = 2. It corresponds to a transfer of 
matter from one density peak towards the other one. The 
end-product of this process will be the equilibrium state 
±1 which is the unique non- linearly stable equilibrium, as 
seen in sect. 14.21 The growth rate of this instability quickly 
goes to zero at low T as the underdense region which con- 
nects both density peaks reaches very low densities from 
ea. (|24f> and the density peaks behave as almost isolated 
subsystems. Besides, it can be seen from the shape of §r) 
in fig. OH that as expected rj + 8r) corresponds to two un- 
equal density peaks at x — 0, L, which follow the profiles 
associated with isolated relaxed density peaks. Thus, in 
agreement with (|23|l we find that in the case displayed 
in fig. 03 the peak at x = which gains some mass gets 
a higher density maximum and a smaller extension while 
the opposite holds for the peak at x = L which looses 
mass. This can be checked more precisely from ea. (l61fl 
as follows. The modified Lane-Emden ea. l|12[) can also be 
written as: 



^2 (-In 77) =7? 



1 , whence 



+ r/=l, 



(62) 



rj_ _ iy_ 
7] rj 2 

which implies that if both 77 and rj + Srj are isothermal 
equilibrium solutions the difference drj satisfies at first or- 
der: 

^-2^'+fl-5 + 2^V^ = 0. (63) 



If' 



On the other hand, taking the derivative of the eigenmode 
equation l|61|) yields (with q' = Srj): 



2 — ,577 - 

77 r/ z 



rf if 
r\ A rf 



Srj 



2II • 64 

77 rj z 



Therefore, for v — we recover ea. (|63|l . In fact, it is 
clear that if (pi,i>i) and (^2,^2) are two equilibrium so- 
lutions of the hydrodynamic eas. H56|l - (|57JI the difference 
(p2 — pi) V2 — Vi) is a static solution at linear order of 
the linearized equations of motion about (pi,v\). Thus, 
in agreement with fig. [5] the unstable eigenmode v\ — > + 
obtained at low temperature for the equilibrium n = 2 
actually yields a quasi-static transfer of matter between 
the two density peaks in the sense that they follow equi- 
librium profiles (parameterized by their varying mass) at 
the fixed temperature T. 

It can be checked by a numerical computation of the 
eigenmode problem H61(l that the behaviours obtained for 
equilibria ±2 extend in a natural manner to higher-order 
equilibrium states ±n. Thus, equilibria ±n show n— 1 un- 
stable eigenmodes (as for the non-linear stability and ther- 
modynamical stability analysis) which can be understood 
as translations of the density peaks or as slow transfers of 
matter between neighbouring density peaks. For instance, 
the equilibrium n — 4 which corresponds to three density 
peaks at x — 0, L/2 and L exhibits one eigenvalue V\ ~ 1 
and two positive eigenvalues > v% close to zero. For 
the mode v\ the central peak rolls to the left or to the 
right. For the mode V2 which is symmetric with respect 
to x — L/2 there is a transfer of matter from the central 
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peak towards both left and right peaks (or the opposite 
since ~q is also a growing eigenmode). For the mode ^3 
which is antisymmetric with respect to x = L/2 (for 8rj) 
there is a transfer of matter from the right peak to the 
central one and a second transfer with the same ampli- 
tude from the central peak to the left one (i.e. a flow from 
the right to the left or the opposite for —q). On the other 
hand, the equilibrium n — —4 which corresponds to two 
density peaks at x = L/4 and 3L/4 exhibits two eigenval- 
ues v\ > vi close to unity and one positive eigenvalue U3 
close to zero. For the mode v\ both peaks move with the 
same velocity towards the same direction while for mode 
V2 they move in opposite directions. The mode 1/3 involves 
a transfer of matter between both peaks. Note that the 
largest growing rate corresponds to both peaks moving 
together in the same direction: the system is not analo- 
gous to two particles in the external potential V alone as 
one needs to take into account the gravitational interac- 
tion $ which balances V at equilibrium. We would obtain 
similar results for higher n. As noticed above, it would be 
interesting to follow the evolution of the system by nu- 
merical simulations to study the possible merging of these 
density peaks. 



5. Vlasov equation 

As discussed in sect. El in the mean field limit where the 
mass of particles goes to zero the iV— body system (J2J is 
actually described by the Vlasov-Poisson system similar to 
0. As is well-known, the equilibrium states obtained in 
the hydrodynamical approach of sect.0]are also stationary 
solutions of the Vlasov-Poisson system with the Maxwell- 
Boltzmann distribution f(x,v) of cq. l|lU|) . Therefore, we 
recover the equilibrium profiles described in sect. I3.1"T1 

On the other hand, as described in Chavanis (2003) 
and Tremaine et al. (1986), the Vlasov dynamics conserves 
any quantity of the form H[f] = — J C(f )dxdv so that 
maxima of any functional H where C(f) is convex are non- 
linearly stable. A particular case of H[f] is the entropy 
S[f] as given in cq.© hence we can already infer that the 
equilibrium ±1 which is a maximum of S is stable below 
T c \. However, this criterion cannot decide whether other 
equilibria ±n with n > 2 are stable or not. 

Next, in order to study the dynamical linear stabil- 
ity of these equilibrium solutions we linearize the equa- 
tions of motion. We first introduce the action-angle vari- 
ables (J, to) which describe the motion of a particle of en- 
ergy E along its orbit in the equilibrium potential (f>o(x) 
from position x- to x+ (see Fridman & Polyachenko 1984; 
Polyachcnko & Shukhman 1981): 



2 r x + . 

J = — / dx^2(E 
2tt J x _ 



<t>a{x)) , w = — 



dS_ 

dJ' 



with: 



S(J,x)= dx'y/2{E-Mx'))- 



(65) 



(66) 



n = — 2 




Fig. 12. The stability eigenvalues aJk/y/SgJ) = iXk/V^9P 
as a function of inverse temperature [3 for the equilibria 
0, ±1 and ±2 for the Vlasov equation of motion. Stability 
corresponds to At < 0. The second eigenvalues A2 ob- 
tained for equilibria ±2 are equal. 

The action-angle variables (J,w) obey Hamilton's equa- 
tions: 

-f7- n - '~fi-* < 67 » 

while the Vlasov equation reads: 
d ^ df df -Of 

■dl=m +w d^ + J dj = - (68) 

We write the distribution function as / = fo(J)+fi(J, w) 
where fo(J) = fo(E) is the equilibrium solution and fi the 
perturbation. Then, the linearized Vlasov equation for fx 
writes: 

d/o 



dfi q dfi_ _ <->h>\ 

dt dw dw d J ' 



(69) 



where 4>\ is the perturbed potential. Because of the period- 
icity of 27T with respect to w of the motion of the particles 
we can write <j)\ as: 



<j)i{j,w,t) = — ^2 iuiJ) 



-i{u)t—lw) 



(70) 



where we look for eigenfrequencies u>. Substituting into 
eq.lJSSJ yields: 



1 



■ilw 



with /r = ^e' 



m 



(71) 



Therefore, Poisson's equation reads: 

lw 

(72) 



0J ■ 



m 



In order to solve the eigenvalue problem l|72|l it is con- 
venient to project this equation onto a biorthonormalized 
basis {p a ,4>a} such that: 

x' = ~3a,a' j (^3) 



= 2g p a and J dx p a t 
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where 5 a ,a' is the usual Kronecker symbol (Fridman & 
Polyachenko 1984; Kalnajs 1977). For instance, we can 
choose with a > 1: 



Pa 



: COS 



^a = -2 



OTT 



which also satisfy the boundary conditions 4>'(0) = 
4>' '(£) = 0. Thus, if we decompose <fii and p\ as <fii = 
J2a a a4>a and pi = J2 a a aPa, substituting into ea. (|72[l 
and integrating over x with the weight <fi a yields: 



to growing modes (instabilities) that vanish for t — > -co. 
For Im(cj) < we need to deform the integration path in 
ea. (|82(l following the usual Landau analysis (e.g. Binney 
& Tremaine 1987; Ichimaru 1973). This leads to an ad- 
cos {air— ^ , (74) ditional contribution due to the residue at x = Qy/2/a. 
^ ' Thus, for Im(w) < the matrix M° Q , is given by the 

sum of the expression l|82|) and of the contribution of the 
residue which reads: 



Im(w) < : Res(M° c 



0a.es' o ^-«we 



2 /a 2 



(83) 



dxdi 
2tt 



l — — oo a' = l 



LU — IQq 



.(75) 



Using the canonical change of variable dxdv — d Jdw and 
next changing from action J to energy E we obtain for an 
isothermal equilibrium distribution / oc e~P E the eigen- 
value problem: 



det(M QQ /(w) - 5 a , a ') 
with: 

[3 



0. 



M aa ,{u) 



3 -JdEf (E)Y: 



to ■ 



(76) 



4a)l(0a')i- (77) 



For the inhomogeneous equilibrium solutions we also need 
to modify the ea. (|78|l in a similar fashion for Im(w) < 0. 
However, since we have no analytical expression for the 
integrand in ea. 178(1 this is more difficult. Thus, for nu- 
merical purposes we simply add the contribution H83|l 
which is correct close to the critical points. Therefore, 
the curves in fig. ^| are not correct for Afc < far 
from the transition T cn but the sign of X k remains ex- 
act. Following Kalnajs (1977), in order to improve the 
convergence with I of the series in ea. l|78[) we can write 
1 2 lu 2 /{1 2 lu 2 - Cu 2 ) = 1 + lu 2 /(1 2 lu 2 - Cu 2 ) which leads to 



M aa t(uj) — A aa i(\) + B aa i, where Cu = i\, 



(84) 



In terms of the dimensionless variables (|14|l we obtain for 
the basis (|?ljl : 



with: 



A , — 

^aa' — 



-4/3 



4/3 



M aa , = / de 



OO 



e e Co 



E /2 - 2 ^ 2 (^MMt.(78) 
1=1 



cl J$ 



de 



OaHVPa.' H 



01,0 



,(85) 



and: 



where ?/>_ is the minimum of V'(C) and with: 
^ dC, cos(a7rC/Ci) cos(lw + ) 



Baa' — 



775" / 7^(0 COS COT — COS Q 7T — 

a« Pel i Cl \ (1/ V Cl 



.(86) 



i>a)l 



Cl 



ay/e — ip 



(79) 



Here we introduced the reduced energy e = f3v 2 /2 + ip of 
the particle, the frequencies Cu and w defined by: 




= , w 




(80) 



and the angular variable w + which describes the orbit in 
the direction of increasing £ from C- up to C+: 



1 



^0 



1 



C- Cl \/e-tp 



TTLUq 



dC 



1 



C_ Cl y/e-i/} 



(81) 



For the homogeneous solution we simply have ip — 0, Cuq = 
y/e, w+ = tt(/(l and (</>„)/ = <5 Ci ,;/(2a^e). Making the 
change of variable x — y/2e we obtain: 



M aa> = *a,a' 



&1 a 2 \/2TrJ- 



dx e~ x2 ' 2 



(82) 



Thus, we can see from eas. <|76[) . (I82II . that the condition 
of marginal stability w = yields again (3 = a 2 f3 c i with 
a = 1,2,.., which agrees with the critical temperatures 
l|17(l obtained in the thermodynamical and hydrodynami- 
cal approaches. The matrix M^ a , as written above is actu- 
ally defined for Im(w) > 0, see eq.{7UJ), which corresponds 



In ea. (|85|l the Kronecker factor ^0 reduces the term I = 
by a factor 2 as compared with 1^0. The faster con- 
vergence as l/l 2 of the expression l|85(l is convenient for 
numerical purposes as it reduces the number of terms to 
be computed. On the other hand, the constant part B aa ', 
which does not depend on A, is computed separately. Note 
that A aa i does not vanish for A = because of the term at 
I = 0. In practice, in order to solve ea. l|76(l we do not need 
to compute the determinant of the matrix M aa ' ■ Starting 
from the transition temperature T c „, where the equilib- 
rium solution ±71 which we are interested in separates from 
the homogeneous equilibrium state, we follow the eigen- 
value Vk of order k of M aa ' and the growth rate Afe of the 
associated eigenmode is set by the condition 2/fe(Afc) = 1 
(e.g. Kalnajs 1977). Moreover, contrary to the usual case 
in plasma physics the eigenfrcquencies lu are imaginary 
so that A is real. This corresponds to purely growing or 
decaying modes without oscillatory behaviours. Thus the 
matrix M aa ' (A) is real and symmetric so that all its eigen- 
values Vk are real, which simplifies the computation of the 
roots of v k (\ k ) = 1. 

We display in fig.^|the stability eigenvalues \ k for the 
equilibria n = 0, ±1 and ±2. We can check that the homo- 
geneous solution n = shows the same stability properties 
as those obtained in the thermodynamical and hydrody- 
namical approaches, as discussed below eq.JHU). On the 
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Fig. 13. The eigenmode fc = 1 for the equilibrium solu- 
tion n — — 2 at temperature (3 = 6/3 c i. We plot both the 
equilibrium profile r\ (dashed line) and its unstable mode 
Srj (solid line) with eigenvalue Ai > 0. 



G 



O 



-2 



-1 













Vlasov 










n = 2, k=l 


/ 




V 




/ 












1 




1 



12 



Fig. 14. The eigenmode fc = 1 for the equilibrium solution 
n = 2 at temperature /3 = 6/3 c i . We plot both the equilib- 
rium profile 77 (dashed line) and its least stable mode 8rj 
(solid line) with eigenvalue Ai ~ 0. 

other hand, from ea. (|82|) we can see that at low temper- 
atures we have for the highest eigenvalue Ai the asymp- 
totic behaviour Ai ~ y/ (3/2f3 c i which yields Ai — + \J2g~p. 
For the equilibrium n — ±1 we only plot the curves up 
to 3/3 c i for clarity and because our formula (|83|l for the 
residue contribution only holds close to T c %. However, we 
checked numerically that A^ remains negative at low T 
(the approximation H83[) only enters for A& < so that 
the sign of the eigenvalues Xk does not depend on it). 
Therefore, as for the thermodynamical and hydrodynam- 
ical approaches the states ±1 are stable below T c \. 

In a similar fashion, the equilibrium n = —2 (which 
corresponds to a density peak at x = L/2) shows one 
unstable mode as in fig. ISJwith \\/ \J2g~p — ► 1 at low tem- 
perature. However, fig. 1131 shows that the corresponding 
eigenmode is significantly different from the one obtained 



in the hydrodynamical framework (compare with fig. I1U|) . 
Indeed, the growing mode br\ no longer tends to the deriva- 
tive 7/ associated with a pure translation of the central 
density peak. Although close to the peak at x = L/2 it can 
be interpreted as a translation the perturbation extends 
as far as the boundaries x = 0, L, so that it also involves a 
deformation of the density profile. Therefore the dynamics 
of the relaxation towards the stable equilibrium ±1 pro- 
ceeds in a manner which is different from the behaviour 
obtained within the hydrodynamical model. In particular, 
the deformation implied by this instability suggests that a 
solitonic behaviour cannot arise as some mass leaks from 
the central peak down to the boundary x — while the 
peak moves towards x = 0. This difference with respect to 
the hydrodynamical case may be understood as follows. 
In the collisional fluid model the evolution proceeds in a 
well-ordered manner through local processes, as each fluid 
element at location x reacts to slow changes in the local 
pressure and density. By contrast, in the collisionless case 
described by the Vlasov equation particles probe far away 
regions as their orbit extends to both sides of the density 
peak. Thus, as the central peak moves towards x — as in 
n g- El the left/right symmetry is broken and the poten- 
tial (j> becomes deeper to the left of the peak. Then, if we 
consider a subsample of particles which relax to the new 
potential as p S ub oc e~X as they can move from one side 
of the density peak to the other, they will populate the 
left side of the peak more than the right side. This yields 
(^77(0) > > 6r)(L) as in fig. ^] (of course the collisionless 
Vlasov dynamics conserves entropy but such a relaxation 
may be applied to the coarse-grained distribution which 
mixes the two subsamples). 

Finally, the equilibrium n = 2 also shows at first an 
unstable mode Ai which quickly goes towards zero as in 
fig. El but contrary to the hydrodynamical case the eigen- 
value Ai becomes negative which means that at low T 
this state turns stable. The density profile of this eigen- 
mode displayed in fig.^Jat 6/3 c i shows that it is again sig- 
nificantly different from its hydrodynamical counterpart. 
Indeed, as for equilibrium —2 it exhibits a broader exten- 
sion than within the hydrodynamical framework which is 
related to the collisionless nature of the system which en- 
tails changes on a more global scale. This implies that 
rj + 5rj does not correspond to a quasi-static evolution 
where both density peaks would follow equilibrium profiles 
with different slowly varying masses. Thus, it appears that 
the decrease of entropy S[f], understood as an H- function 
as in Chavanis (2003), due to the distortion of the pro- 
files close to the boundaries x = 0, L, farther from iso- 
lated equilibrium states dominates the increase of entropy 
which would be associated with the transformation of the 
system from two separate density peaks towards only one 
density peak (state ±1 which exhibits a higher S[f], as 
seen from fig. |2J). Thus, thermodynamical considerations 
or non-linear analysis based on H-functions are not suffi- 
cient to obtain the dynamical stability properties of the 
system. One needs to take into account the constraints 
brought by the Vlasov dynamics: one cannot select any 
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arbitrary deformation path from equilibrium 2 towards 
equilibrium ±1. As discussed in Lynden-Bell (1967), the 
infinite number of conserved quantities associated with 
the collisionless dynamics (the total mass of phase-space 
element greater than any / is conserved) can actually be 
taken into account by using a new "collisionless entropy" 
In W which differs from the usual Boltzmann entropy. On 
the other hand, the comparison with sect. 0] shows that 
the collisional and collisionless dynamics are significantly 
different, as it appears that the hydrodynamical dynam- 
ics allows the system to flow from state n — 2 to state 
n = ±1 (which implies that equilibrium n — 2 is unsta- 
ble) whereas at low T the Vlasov dynamics prevents the 
system from going to n = ±1 from n — 2 through lin- 
ear instabilities (but clearly a non-linear perturbation can 
lead to relaxation towards n = ±1). 

6. Conclusion 

We have studied in this article the thermodynamical prop- 
erties of a 1-D gravitational system built from the prob- 
lem of large-scale structure formation in cosmology. We 
have considered time-scales much smaller than the Hubble 
time so that the expansion of the universe can be ne- 
glected and the system can be restricted to a finite size 
L. Then, the cosmological background merely yields an 
effective external potential V in addition to the gravita- 
tional self-interaction $ so that an homogeneous equilib- 
rium (corresponding to the Hubble flow) exists at all tem- 
perature T. We have found that all three micro-canonical, 
canonical and grand-canonical ensembles give the same 
results. There exists a series of critical temperatures T cn 
(T c i > T C 2 > ..) below which two new equilibria ±n 
appear, reflecting the usual Jeans instability mechanism. 
Above T c i the homogeneous state is thermodynamically 
stable whereas below T c \ it is unstable and only equilib- 
ria ±1 (corresponding to a density peak at the boundary 
x — or L) are stable while other states ±n with n > 2 are 
unstable. The system shows a second-order phase transi- 
tion at T c i where the specific heat is discontinuous. 

Next, we have studied the hydrodynamical properties 
of these equilibria. We have found that both the non-linear 
stability and the dynamical linear stability coincide with 
the results of the thermodynamical analysis. The evolution 
of the system can be understood in a simple manner from 
the translation and diffusion of the various density peaks 
associated with states ±n which would relax towards the 
stable equilibria ±1. Finally, we have studied the colli- 
sionless dynamics governed by the Vlasov equation. The 
equilibria ±1 are again both linearly and non-linearly sta- 
ble while the homogeneous state is unstable below T cl . We 
found that the equilibrium n = — 2 (one density peak at 
x — L/2) is again unstable but the equilibrium n — 2 (two 
density peaks located at both boundaries) shows a new be- 
haviour. Although it is unstable just below T C 2 as in the 
thermodynamical and hydrodynamical approaches it now 
turns stable at low temperatures. This departure from the 
results of the thermodynamical and hydrodynamical anal- 



ysis is due to the dynamical constraints associated with 
the equations of motion of a collisionless system. On the 
other hand, other equilibria ±n with n > 3 are again un- 
stable. Therefore, depending on the initial conditions the 
system would relax towards equilibria ±1 or 2 (i.e. a den- 
sity peak at either one or both boundaries). 

Thus, the behaviour of the 1-D gravitational system 
(OSC) studied in this article is quite similar to the HMF 
model, described for instance in Chavanis et al. (2005). 
We obtain the same second-order phase transition from 
a high-temperature homogeneous equilibrium to a clus- 
tered equilibrium ±1 at low temperature. The stability 
properties of these two configurations are also similar, for 
all three thermodynamical, hydrodynamical and collision- 
less analysis. However, in the case of the OSC system the 
scale-free nature of the gravitational interaction leads to 
an additional series of equilibria ±n (n > 2). These new 
states are unstable, except for the equilibrium n = 2 which 
turns linearly stable for a collisionless dynamics at low T, 
but remains unstable for both the thermodynamical and 
hydrodynamical approaches. This provides an interesting 
case where the collisionless dynamics exhibits a specific 
behaviour and leads to a relaxation which clearly depends 
on initial conditions. On the other hand, it would be in- 
teresting to study the non-linear evolution of the system 
and the possible merging of several density peaks in both 
the hydrodynamical and collisionless cases. This might be 
studied through N— body simulations as in Hohl & Feix 
(1967) or through a direct computation of the Vlasov dy- 
namics as in Alard & Colombi (2005). Finally, one should 
include the expansion of the universe to study the non- 
equilibrium dynamics associated with cosmological hier- 
archical scenarios. This is left for future studies but we 
can expect this 1-D gravitational system to provide a use- 
ful tool to investigate such processes. 
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